real function integral1(f,a,b)
    real h
    integer i,n
    n=1000
    h=(b-a)/real(n)
    integral1=0.5*(f(a)+f(b))
    do, i=1,n-1
        integral1=integral1+f(a+i*h)
    end do
    integral1=h*integral1
end function

real function integral2(f,a,b)
    real h,f1,f2
    integer i,n
    n=2
    f1=0
    f2=4!f2-f1>1e-5
    do while(abs(f1-f2) .GT. 1e-5)
        f2=f1
        n=n*2
        f1=(f(a)+f(b))/2.0
        h=(b-a)/real(n)
        do i=1,n-1
            f1=f1+f(a+i*h)
            !print *,'time=',i
        end do
        f1=f1*h
        !print*, 'f1=',f1
    end do
    integral2=f1
end function

real function y1(x)
    y1=sqrt(1+sin(x))
end function

real function y2(x)
    y2=x*x+x+1
end function

real function y3(x)
    y3=exp(-0.5*x**2)
end function

program chapter11_13
    implicit none
    real integral1, integral2, y1, y2, y3
    external y1,y2,y3
    intrinsic sin,sqrt,exp
    print* , 'The first method'
    print *, integral1(y1,0.0,3.1415926/2.0)
    print *, integral1(y2,0.0,2.0)
    print *, integral1(y3,0.0,1.0)
    print *, 'The second method'
    print *, integral2(y1,0.0,3.1415926/2.0)
    print *, integral2(y2,0.0,2.0)
    print *, integral2(y3,0.0,1.0)
end program
